function   model_sol = find_solution( p_2D_init,  par_input  )
% Usage: find_solution( p_2D_init,  par  ).
% 
% Notes: This function solves the model with parameters alpha, beta, and initial values p_2D_init. p_2D_init is a 2D matrix.
% All other parameters are stored in par.   
% This function will call "find_solution_inner_iteration" to generate a new price function from an old price function. 

% status = evalWithTimer( 'par = initialization(par_input)', 3);
par = initialization(par_input);

%% *** Solve the model ***
tol = par.tol;
errors = [];  iter=0;  max_iter=15;
p_2D_new = p_2D_init;    dampen = 0.5 ;  % Note: dampening the iteration is very important. It improves both stability and convergence.
tic
while(iter<max_iter)
    disp(  [  '************ Solve the 2D problem, iter=', num2str(iter),  ' for ', par.model_name, ' model ************' ]    )
    p_2D_last = p_2D_new;
    model_sol =  find_solution_inner_iteration( p_2D_new,  par  );
    p_2D_new = (1-dampen)*model_sol.p_matrix + dampen*p_2D_last;
    
    iter = iter + 1;
    errors(iter) = sum(sum( abs(  p_2D_new - p_2D_last ) ));
    disp( [  '    error =  ', num2str(errors(iter)) ] )
    if( errors(iter) < tol )
        break;
    end
    if( length(errors)>=2 && log(errors(end)/errors(end-1))>0.5 )   % if there is coniderable amount of reversion in error.
        break;
    end
    toc
end


%% Resolve kappap
kappap_matrix =  resolving_kappap( p_2D_new,  par  );


%% Output
model_sol.kappap_matrix = kappap_matrix;
model_sol.par = par;
model_sol.errors = errors;

